/**
 * @file core.h
 * 
 * @brief ...
 */
#ifndef _CORE_H_
#define _CORE_H_

#include <functional>
#include <algorithm>
#include <iostream>
#include <utility>
#include <cassert>

#define NAMESPACE_LA_BEGIN namespace RubbishLinearAlgebra {
#define NAMESPACE_LA_END }

NAMESPACE_LA_BEGIN

/**
 * @class SharedView
 * 
 * @brief ...? 
 */
template <typename T>
struct SharedView {
    int N;
    T *data;

    SharedView(int n, T* ptr) : N(n), data(ptr) {}
};

/**
 * @class Matrix
 * 
 * @brief a naive row-major m by n matrix
 */
template <typename T>
struct Matrix {
private:
    int M, N;
    T *data;

public:
    Matrix() : M(0), N(0), data(nullptr) {}

    Matrix(int m, int n) : M(m), N(n)
    {
        assert(m > 0 && n > 0);
        data = new T[m * n];
    }

    Matrix(int m, int n, const SharedView<T>& view) : M(m), N(n)
    {
        assert(m > 0 && n > 0);
        assert(m * n == view.N);
        data = view.data;
    }

    Matrix(const Matrix& rhs)
    {
        M = rhs.M, N = rhs.N;
        if (rhs.data != nullptr) {
            data = new T[M * N];
            std::copy(rhs.data, rhs.data + M * N, data);
        } else {
            data = nullptr;
        }
    }

    Matrix(Matrix&& rhs)
    {
        M = rhs.M, N = rhs.N;
        data = rhs.data;
        rhs.M = rhs.N = 0;
        rhs.data = nullptr;
    }

    Matrix& operator = (const Matrix& rhs)
    {
        M = rhs.M, N = rhs.N;
        if (rhs.data != nullptr) {
            data = new T[M * N];
            std::copy(rhs.data, rhs.data + M * N, data);
        } else {
            data = nullptr;
        }
        return *this;
    }

    Matrix& operator = (Matrix&& rhs)
    {
        M = rhs.M, N = rhs.N;
        data = rhs.data;
        if (this != &rhs) {
            rhs.M = rhs.N = 0;
            rhs.data = nullptr;
        }
        return *this;
    }

    ~Matrix()
    {
        if (data != nullptr) {
            delete[] data;
            data = nullptr;
        }
    }

    int m() const { return M; }
    int n() const { return N; }

    std::pair<size_t, size_t> size() const
    {
        return std::make_pair<size_t, size_t>(M, N);
    }

    SharedView<T> view() const
    {
        return SharedView<T>(M * N, data);
    }

    /**
     * @warning destroys all elements!
     */
    void resize(int m, int n)
    {
        assert(m > 0 && n > 0);
        if (data != nullptr)
            delete[] data;
        M = m, N = n;
        data = new T[M * N];
    }

    T& operator () (int i, int j)
    {
        assert(0 <= i && i < M && 0 <= j && j < N);
        assert(data != nullptr);

        return data[i * N + j];
    }

    const T& operator () (int i, int j) const
    {
        assert(0 <= i && i < M && 0 <= j && j < N);
        assert(data != nullptr);

        return data[i * N + j];
    }

    friend std::ostream& operator << (std::ostream& os, const Matrix<T>& mat)
    {
        os << "matrix ( " << mat.m() << " x " << mat.n() << " ):\n";
        for (int i = 0; i < mat.m(); ++i) {
            os << '\t';
            for (int j = 0; j < mat.n(); ++j) {
                os << mat(i, j) << ' ';
            }
            os << std::endl;
        }
        return os;
    }

    /** @note: inplace operations **/
    void zero();
    void eye(const T& v);
    void identity();
    // void transpose();

    Matrix& operator += (const T& v);
    Matrix& operator -= (const T& v);
    Matrix& operator *= (const T& v);
    Matrix& operator /= (const T& v);
    Matrix& operator += (const Matrix& rhs);
    Matrix& operator -= (const Matrix& rhs);
    Matrix& operator *= (const Matrix& rhs);
    Matrix& operator /= (const Matrix& rhs);

    Matrix operator - () const;
    Matrix operator + (const T& v) const;
    Matrix operator - (const T& v) const;
    Matrix operator * (const T& v) const;
    Matrix operator / (const T& v) const;
    Matrix operator + (const Matrix& rhs) const;
    Matrix operator - (const Matrix& rhs) const;
    Matrix operator * (const Matrix& rhs) const;
    Matrix operator / (const Matrix& rhs) const;

    void from_istream(std::istream& is);
};

/**
 * @class Vector
 * 
 * @brief a naive Vector
 */
template <typename T>
struct Vector {
private:
    int N;
    T* data;

public:
    Vector() : N(0), data(nullptr) {}
    
    Vector(int n) : N(n)
    {
        assert(n > 0);
        data = new T[n];
    }

    Vector(int n, const SharedView<T>& view) : N(n)
    {
        assert(n > 0);
        assert(n == view.N);
        data = view.data;
    }

    ~Vector()
    {
        if (data != nullptr) {
            delete[] data;
            data = nullptr;
        }
    }

    int n() const { return N; }
    size_t size() const { return N; }

    SharedView<T> view() const
    {
        return SharedView<T>(N, data);
    }

    T& operator [] (int i)
    {
        assert(0 <= i && i < N);
        assert(data != nullptr);

        return data[i];
    }

    const T& operator [] (int i) const
    {
        assert(0 <= i && i < N);
        assert(data != nullptr);

        return data[i];
    }

    friend std::ostream& operator << (std::ostream& os, const Vector<T>& vec)
    {
        os << "vector ( " << vec.n() << " ):\n";
        os << '\t';
        for (int i = 0; i < vec.n(); ++i)
            os << vec[i] << ' ';
        os << std::endl;
        return os;
    }

    Vector& operator += (const T& v);
    Vector& operator -= (const T& v);
    Vector& operator *= (const T& v);
    Vector& operator /= (const T& v);
    Vector& operator += (const Vector& rhs);
    Vector& operator -= (const Vector& rhs);
    Vector& operator *= (const Vector& rhs);
    Vector& operator /= (const Vector& rhs);

    Vector operator - () const;
    Vector operator + (const T& v) const;
    Vector operator - (const T& v) const;
    Vector operator * (const T& v) const;
    Vector operator / (const T& v) const;
    Vector operator + (const Vector& rhs) const;
    Vector operator - (const Vector& rhs) const;
    Vector operator * (const Vector& rhs) const;
    Vector operator / (const Vector& rhs) const;

    void from_istream(std::istream& is);
};

// matrix arithmetic
template <typename T>
Matrix<T> Matrix<T>::operator - () const
{
    assert(data != nullptr);
    Matrix<T> res(M, N);
    std::transform(data, data + M * N, res.data, std::negate<T>());
    return res;
}

template <typename T>
Matrix<T>& Matrix<T>::operator += (const T& v)
{
    assert(data != nullptr);
    std::transform(data, data + M * N, data, [&v](const T& u) -> T { return u + v; });
    return *this;
}

template <typename T>
Matrix<T>& Matrix<T>::operator -= (const T& v)
{
    assert(data != nullptr);
    std::transform(data, data + M * N, data, [&v](const T& u) -> T { return u - v; });
    return *this;
}

template <typename T>
Matrix<T>& Matrix<T>::operator *= (const T& v)
{
    assert(data != nullptr);
    std::transform(data, data + M * N, data, [&v](const T& u) -> T { return u * v; });
    return *this;
}

template <typename T>
Matrix<T>& Matrix<T>::operator /= (const T& v)
{
    assert(data != nullptr);
    std::transform(data, data + M * N, data, [&v](const T& u) -> T { return u / v; });
    return *this;
}

template <typename T>
Matrix<T> Matrix<T>::operator + (const T& v) const
{
    assert(data != nullptr);
    Matrix<T> res(M, N);
    std::transform(data, data + M * N, res.data, [&v](const T& u) -> T { return u + v; });
    return res;
}

template <typename T>
Matrix<T> Matrix<T>::operator - (const T& v) const
{
    assert(data != nullptr);
    Matrix<T> res(M, N);
    std::transform(data, data + M * N, res.data, [&v](const T& u) -> T { return u - v; });
    return res;
}

template <typename T>
Matrix<T> Matrix<T>::operator * (const T& v) const
{
    assert(data != nullptr);
    Matrix<T> res(M, N);
    std::transform(data, data + M * N, res.data, [&v](const T& u) -> T { return u * v; });
    return res;
}

template <typename T>
Matrix<T> Matrix<T>::operator / (const T& v) const
{
    assert(data != nullptr);
    Matrix<T> res(M, N);
    std::transform(data, data + M * N, res.data, [&v](const T& u) -> T { return u / v; });
    return res;
}

template <typename T>
Matrix<T> Matrix<T>::operator + (const Matrix<T>& rhs) const
{
    assert(M == rhs.M && N == rhs.N && data != nullptr);
    Matrix<T> res(M, N);
    std::transform(data, data + M * N, rhs.data, res.data, std::plus<T>());
    return res;
}

template <typename T>
Matrix<T> Matrix<T>::operator - (const Matrix<T>& rhs) const
{
    assert(M == rhs.M && N == rhs.N && data != nullptr);
    Matrix<T> res(M, N);
    std::transform(data, data + M * N, rhs.data, res.data, std::minus<T>());
    return res;
}

template <typename T>
Matrix<T> Matrix<T>::operator * (const Matrix<T>& rhs) const
{
    assert(M == rhs.M && N == rhs.N && data != nullptr);
    Matrix<T> res(M, N);
    std::transform(data, data + M * N, rhs.data, res.data, std::multiplies<T>());
    return res;
}

template <typename T>
Matrix<T> Matrix<T>::operator / (const Matrix<T>& rhs) const
{
    assert(M == rhs.M && N == rhs.N && data != nullptr);
    Matrix<T> res(M, N);
    std::transform(data, data + M * N, rhs.data, res.data, std::divides<T>());
    return res;
}

// vector arithmetic
template <typename T>
Vector<T> Vector<T>::operator - () const
{
    assert(data != nullptr);
    Vector<T> res(N);
    std::transform(data, data + N, res.data, std::negate<T>());
    return res;
}

template <typename T>
Vector<T>& Vector<T>::operator += (const T& v)
{
    assert(data != nullptr);
    std::transform(data, data + N, data, [&v](const T& u) -> T { return u + v; });
    return *this;
}

template <typename T>
Vector<T>& Vector<T>::operator -= (const T& v)
{
    assert(data != nullptr);
    std::transform(data, data + N, data, [&v](const T& u) -> T { return u - v; });
    return *this;
}

template <typename T>
Vector<T>& Vector<T>::operator *= (const T& v)
{
    assert(data != nullptr);
    std::transform(data, data + N, data, [&v](const T& u) -> T { return u * v; });
    return *this;
}

template <typename T>
Vector<T>& Vector<T>::operator /= (const T& v)
{
    assert(data != nullptr);
    std::transform(data, data + N, data, [&v](const T& u) -> T { return u / v; });
    return *this;
}

template <typename T>
Vector<T> Vector<T>::operator + (const T& v) const
{
    assert(data != nullptr);
    Vector<T> res(N);
    std::transform(data, data + N, res.data, [&v](const T& u) -> T { return u + v; });
    return res;
}

template <typename T>
Vector<T> Vector<T>::operator - (const T& v) const
{
    assert(data != nullptr);
    Vector<T> res(N);
    std::transform(data, data + N, res.data, [&v](const T& u) -> T { return u - v; });
    return res;
}

template <typename T>
Vector<T> Vector<T>::operator * (const T& v) const
{
    assert(data != nullptr);
    Vector<T> res(N);
    std::transform(data, data + N, res.data, [&v](const T& u) -> T { return u * v; });
    return res;
}

template <typename T>
Vector<T> Vector<T>::operator / (const T& v) const
{
    assert(data != nullptr);
    Vector<T> res(N);
    std::transform(data, data + N, res.data, [&v](const T& u) -> T { return u / v; });
    return res;
}

template <typename T>
Vector<T> Vector<T>::operator + (const Vector<T>& rhs) const
{
    assert(N == rhs.N && data != nullptr);
    Vector<T> res(N);
    std::transform(data, data + N, rhs.data, res.data, std::plus<T>());
    return res;
}

template <typename T>
Vector<T> Vector<T>::operator - (const Vector<T>& rhs) const
{
    assert(N == rhs.N && data != nullptr);
    Vector<T> res(N);
    std::transform(data, data + N, rhs.data, res.data, std::minus<T>());
    return res;
}

template <typename T>
Vector<T> Vector<T>::operator * (const Vector<T>& rhs) const
{
    assert(N == rhs.N && data != nullptr);
    Vector<T> res(N);
    std::transform(data, data + N, rhs.data, res.data, std::multiplies<T>());
    return res;
}

template <typename T>
Vector<T> Vector<T>::operator / (const Vector<T>& rhs) const
{
    assert(N == rhs.N && data != nullptr);
    Vector<T> res(N);
    std::transform(data, data + N, rhs.data, res.data, std::divides<T>());
    return res;
}


// matrix multiplication
template <typename T>
Matrix<T> matmul(const Matrix<T>& lhs, const Matrix<T>& rhs)
{
    assert(lhs.n() == rhs.m());
    int M = lhs.m(), N = rhs.n(), P = lhs.n();
    Matrix<T> res(M, N);
    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < N; ++j) {
            T v = 0;
            for (int k = 0; k < P; ++k)
                v += lhs(i, k) * rhs(k, j);
            res(i, j) = std::move(v);
        }
    }
    return res;
}

template <typename T>
T dot(const Vector<T>& lhs, const Matrix<T>& rhs)
{
    assert(lhs.n() == rhs.n());
    T res = 0;
    for (int i = 0; i < lhs.n(); ++i)
        res += lhs[i] * rhs[i];
    return res;
}

// matrix initialization
template <typename T>
void Matrix<T>::zero()
{
    assert(M > 0 && N > 0 && data != nullptr);
    std::fill(data, data + M * N, T(0));
}

template <typename T>
void Matrix<T>::eye(const T& v)
{
    zero();
    for (int i = 0; i < std::min(M, N); ++i)
        data[i * N + i] = v;
}

template <typename T>
void Matrix<T>::identity()
{
    eye(1);
}

// file I/O
template <typename T>
void Matrix<T>::from_istream(std::istream& is)
{
    assert(M != 0 && N != 0);
    for (int i = 0; i < M * N; ++i)
        is >> data[i];
}

template <typename T>
void Vector<T>::from_istream(std::istream& is)
{
    assert(N != 0);
    for (int i = 0; i < N; ++i)
        is >> data[i];
}

NAMESPACE_LA_END

#endif // _CORE_H_